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Abstract. 

Upgraded spectroscopic hardware and an improved impurity concentration calculation allow accurate 
determination of boron density in the ASDEX Upgrade tokamak. A database of boron measurements is 
compared to quasilinear and nonlinear gyrokinetic simulations over a range of H-mode plasma regimes. 
. The peaking of the measured boron profiles shows a strong anti-correlation with the plasma rotation, via a 

■ relationship explained and reproduced by the theory. It is demonstrated that the rotodiffusive impurity flux 

driven by the rotation gradient is required for the modelling to reproduce the hollow boron profiles at higher 
rotation, and that each symmetry breaking mechanism that causes momentum transport also couples to 
rotodiffusion. At lower rotation, the parallel compressive convection is required to match the most peaked 
boron profiles. The sensitivities of both datasets to possible errors is investigated, and quantitative agreement 
* ^ 1 is found within the estimated uncertainties. The approach used can be considered a template for mitigating 

uncertainty in quantitative comparisons between simulation and experiment. 



PACS numbers: 52.25.Fi, 52. 55. Fa, 52.65.Tt, 52.25.Vy, 52.35.Ra 



1. Introduction 

To achieve optimum fusion performance, future reactors will need to control the build up of impurities in 
the plasma core. Light impurities dilute the fusion fuel, while heavier impurities determine radiation losses 
and impact plasma control. The design and operation of future reactors will therefore benefit greatly from 
predictive capability for modelling impurity transport under different plasma conditions. 

Core impurity transport is observed to be anomalous in both H modes and L modes [TJ [5], though 
neoclassical transport can also contribute, particularly for heavier impurities and at higher collisionalities. 
In the edge, transport of impurities has been observed to follow the neoclassical predictions [3]. The tools 
for modelling both neoclassical and turbulent impurity transport have reached the stage of development 
at which it is possible to make comparisons with experiment. Quantitative comparisons across a range of 
plasma conditions are necessary to validate the appropriateness of these tools before they can be applied 
with confidence to make predictions for future devices [1]. 

Turbulent transport is well known to be stiff with respect to the drive, making quantitative comparison 
of fluxes challenging in all transport channels. This difficulty can be avoided by the examination of steady 
state phases in which dimensionless ratios of transport coefficients are used to predict null flux density 
gradients which can be compared directly with the measurements [5]. For impurities, one obstacle to such 
efforts is the difficulty of making accurate measurements of impurity density. 

The aim of this work is to quantitatively compare the best available measurements to physically 
comprehensive gyrokinetic and neoclassical modelling. The gyrokinetic modelling is performed with the 
GKW code [6] including the rotodiffusive contribution [7| (the off-diagonal coupling between rotation and 
particle transport channels) and centrifugal effects [8] . The neoclassical transport is computed with the NEO 
code Uni HJ (which also contains the effects of rotation) , but is found to be negligible for the database 
considered. We exploit the recently upgraded charge exchange (CX) recombination spectroscopy diagnostic 
on the ASDEX Upgrade tokamak (AUG) combined with recent improvements in the density interpretation 
|12) to build a database of boron profiles for H-mode plasmas under a variety of heating regimes. This work 
builds on the work of Ref. [5] by adding a more diverse database covering a larger parameter range, with 
more accurate measurements. The present modelling furthermore includes centrifugal effects and a detailed 
investigation of the sensitivities of the results to enable quantitative conclusions. We show that, within the 
uncertainties, the primary correlation in the database, between plasma rotation and impurity peaking, is 
quantitatively reproduced and explained by the modelling. 

This paper is structured as follows: Section II describes the improvements to the CX measurement 
of boron density, Section III presents the experimental database, and the investigation of systematic 
measurement uncertainties. Section IV describes the modelling, including sensitivities studies, to enable 
a robust quantitative comparison. Conclusions are drawn in Section V. 

2. Measurement of Boron Density 

The core charge exchange (CX) diagnostic on AUG was substantially upgraded prior to the 2011 campaign to 
a 30 channel system which provides much higher temporal resolution and more complete profiles as compared 
with the previous system |13) |14) . This system is routinely used to measure intrinsic boron (the tungsten wall 
is regularly boronised) to provide high resolution measurements of ion temperature and rotation. Impurity 
densities can also be calculated from the CX spectra, but this analysis is more challenging, not only because 
spectral intensities must be correctly calibrated, but also because the subsequent calculation requires correct 
treatment of the neutral beam attenuation and halo contributions up to at least the second excited state [T2] . 
The CHICA Charge exchange Impurity Concentration Analysis code used for these calculations has recently 
been substantially rewritten to model these contributions more accurately, and now includes an improved 
beam geometry based on beam emission spectroscopy measurements, more accurate modelling of the beam 
excited state populations, a calculation of the halo population through the Monte-Carlo fast ion simulation 
code (f90)fidasim jT5J [TB] , and inclusion of the halo contribution to the CX signal using ADAS (Atomic 
Data and Analysis Structure) thermal charge exchange cross-sections [17l [18] . The overall result of these 
improvements has been to reduce the calculated boron densities to ~ 40% of their value in the previous 
simpler analysis. Throughout this paper, when we refer to "boron density measurement" to distinguish it 
from the independent gyrokinetic simulation, the reader should bear in mind that the measurement technique 
requires absolute calibration, multiple sources of atomic data, and a Monte-Carlo simulation of the neutral 
beam. 
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The phases and profiles in a typical shot are shown in Fig. [TJ which provides an example of a somewhat 
generic result for AUG: When electron heating is applied to an NBI heated H-mode, increased density 
peaking and decreased rotation are observed |19) . 

3. Experimental database 

For this work, a database of boron density profiles from H-mode shots from the 2011 campaign on AUG 
has been constructed. The database consists of 29 steady phases of at least 0.25 seconds from 10 shots; 
the energy confinement time for these phases ranges from 0.06 to 0.12 seconds. All shots have a continuous 
2.5 MW NBI heating power from the beam used for the CX measurements, which is predominantly radial 
(major radius of tangency / major radius of axis = Rt / Ra — 0.31). This heating provides a minimum to 
which additional NBI and ECRH heating are added (shots with ICRH are excluded). The database covers a 
range of plasmas with plasma currents I p from 0.6 to 1.0 MA, NBI powers between 2.5 and 7.5 MW, ECRH 
powers from to 2.5 MW, core electron densities between 5 and 12 •10 _19 m -3 , Greenwald fraction n e /naw 
= 0.35-0.85, /3n = 0.65-2.0, and safety factor (795 = 3.9-7.1. All plasmas are lower single null with toroidal 
field -2.5T and similar edge elongation (1.6) and average triangularity (0.15). This database covers a greater 
range of NBI heating and Mach numbers than Ref. [5] which also used earlier charge exchange hardware 
and did not use the improved analysis for the impurity concentration. 

A number of strong correlations are immediately visible in the database (Fig. [2]). The most striking is 
between the Mach number u = v tor /vth and the logarithmic gradient of the boron density, R/L u b (Fig 2g). 
The Mach number is defined with the species thermal velocity Vth — y/2T/m. A similar correlation can also 
be seen between R/L n B and v! = — (R/vth^du/dr since this quantity is highly correlated with the Mach 
number (Fig 2h, where u' B is normalized with the boron thermal velocity). Since temperature gradients 
are the primary quantity in determining the characteristics of turbulent modes, the correlations between 
temperature gradients and R/L n s provide a strong signature that turbulence is important in the impurity 
transport. Clearly, many of the "independent variables" plotted on the x-axes also have correlations between 
them (e.g. heating powers, temperature gradients and rotation): A strong correlation between rotation and 
ion temperature gradient (not shown), is present (since the AUG NBI heating systems always drive a co- 
current rotation), but contains two branches due to different injection angles and deposition radii of the 
NBI beams. The different geometries of the NBI beams are also responsible for the two branches visible in 
the plots of R/L nB against R/L Ti and T e /Tj and P^c/Ptat (Figs, 2f, 2b, 2e, 2c respectively); the points 
with more than 2.5 MW NBI in the highlighted branch all include the most tangential beam (which has 
P-t/Ra = 0.75). Fig. 2d indicates that the safety factor is not a primary dependence, which suggests that 
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Figure 2. Correlations in the AUG experimental H-mode database of boron density measurements at 
r/a = 0.5. The color of the points indicates the ECRH heating power between and 2.5 MW, the shape of 
the points represents the NBI heating power between 2.5 and 7.5 MW. For the points inside or below the 
dashed line, the most tangential NBI beam provides NBI heating in addition to the base 2.5 MW. 



neoclassical transport is not dominant. 

The primary correlation, between rotation and impurity peaking is not caused by a single mechanism, 
but is a non-trivial combination of the correlations between heating powers, temperature gradients, and 
rotation (some of which are a result of the beam geometries specific to AUG). These quantities together 
determine the turbulence characteristics (as will be demonstrated by the modelling in the following sections) 
in such a way that the overall result produces the strong correlation shown. 

3.1. Uncertainties in the measurement 

The comparison of the simulations with experiment over a sufficiently large database helps to avoid random 
errors in the measurements from distorting the global comparison over the database. Systematic errors, by 
contrast, are much harder both to estimate and to eliminate, and can certainly affect the quantitative 
comparison between measurement and modelling. In this section we identify and estimate the largest 
sources of known uncertainty in the measurement; the sensitivities of the modelling to input uncertainties is 
investigated in Sec. 14.11 

The sample errors bars shown in Fig. [T]are the mean error bars provided by the CHICA code and include 
the CX fitting error, in addition to systematic uncertainty estimates for the atomic data and intensity 
calibration. One further uncertainty in the calculation is the deuteron density profile required for the 
fidasim modelling of the neutral contributions, a quantity not routinely measured on AUG. The comparison 
of the FIDASIM modelling with beam emission spectroscopy in Ref. 12 demonstrated that the modelling of 
the halo contribution to the CX signal is accurate to within 30% , part of the discrepancy could be due to 
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Figure 3. Estimated uncertainties in measured R/L„b at mid radius for the entire database. (Left) 
Boron density gradient sensitivity via the halo contribution to different deuterium concentrations. (Right) 
Estimated error bars on the reference points (83% D) combining both systematic and random uncertainty 
estimates using a distribution of spline fits. 



uncertainties in the deuteron density. 

The deuteron profile must be estimated from the electron density profile combined with an assumption 
about the overall impurity concentration. In addition, the profiles of electron density (primarily derived from 
an inversion of 5 horizontal line integrated interferometry measurements) also contribute to the uncertainty. 
The electron profiles used in this database are obtained from an inversion of the interferometry including the 
edge lithium beam, obtained from the Bayesian integrated data analysis (IDA) tool [201 HI] ■ The deuteron 
density influences the beam penetration and attenuation, and the relationship between the deuteron profile 
uncertainty and the calculated boron density cannot be surmised a priori. 

To quantify the influence of the deuteron density uncertainty, three points from the database were 
investigated by calculating the boron density with a range of possible deuterium concentrations. The 
deuterium concentration was varied between {91% , 83% and 78% }, with the latter corresponding to a 
very dirty plasma for AUG (with Z cff = 2.4 assuming 4% He, 1% B, 1% C and 0.04% W). The effect of these 
deuterium concentrations on the calculated boron density gradient is shown for the entire database in Fig. 
|3l The results show that the use of the higher (lower) deuterium concentrations systematically increases 
(decreases) the measured R/L n B at mid radius, due to the variation of the halo contribution to the CX 
signal. We also investigated different inversions for the n e profile and the effect of non-flat Z e g profile; the 
results do not differ significantly from the set of three deuterium concentrations presented. 

To calculate an uncertainty on the logarithmic gradient R/L n s = —{R/nB){driB / dr) we separate 
the global systematic and local random errors for appropriate uncertainty propagation into the gradient 
quantities. Since each measurement is a steady state phase of between 0.25 and 1.0 seconds, the fitting error 
estimate est is the standard error (~ 1/ yjnt) of the CX fitting errors for each time point (where nt between 
3 and 15 is the number of measurements in the time window). The standard deviation of the measurements 
in the time window is considered to be a fluctuation uncertainty £3 and we furthermore define a per radial 
point calibration error £p pce which is the distance from the fitted trend line (up to a maximum of 10% of 
the density value). For each radial point, these uncertainties are combined as 



: fit + £ ppco + e fl' (1) 

Using these radial point errors to add Gaussian noise to each radial point, 10,000 cubic spline fits are 
generated for every point in the database. The distribution of R/L n B from the spline fits is itself Gaussian, 
with mean and median identical to the spline fits on the raw data. The standard deviation of the logarithmic 
gradient of the fits e r i n is combined with the sensitivity to the deuteron profile s nD to produce the error bars 



£tot = Y e rin ~^ £ nn snown in Fig. [3] (global systematic uncertainties, e.g. in global calibration or atomic 
data do not contribute to the logarithmic gradient). 

The size of the gradient uncertainties is influenced by the tension on the spline fits. The value used 
was chosen (prior to the generation of the randomised fits) as the value which best preserved the trend of 
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the data while removing the influence of the features of individual points. While the reader may note that 
there is some subjectivity to this choice, the above exercise demonstrates that fitted logarithmic gradients 
can be much more uncertain than the absolute density at a given location. However, the advantage of 
this dimensionless quantity is that it can be compared with quasi-linear simulations at a single radius (as 
discussed in the next section). This allows a quantitative comparison to be made without a calculation of 
the nonlinear saturation of stiff turbulent fluxes, which are extremely sensitive to the input gradients. 

4. Quasi-linear Gyrokinetic modelling, with neoclassics 

To validate our understanding of the physical mechanism responsible for the trends in the database, 
quasilinear simulations were conducted for each point in the database with the flux tube version of the 
gyrokinetic code GKW [BJ. The boron is modelled in the trace limit in which the microinstabilities are 
completely determined by the bulk plasma conditions. In this limit, the response of impurities to the 
turbulence is completely linear in the gradients of the impurities, even for nonlinear simulations [221 251 124). 
The impurity flux can then be exactly decomposed as 



where the individual terms correspond to the diffusive, thermodiffusive, rotodiffusive, and convective parts, 
respectively. Here we define the rotodiffusive coefficient C u using the impurity rotation gradient v! s 
normalized with the velocity of the bulk ions. The rotodiffusive term [7] is a result of symmetry breaking 
mechanisms and is therefore closely connected with the mechanisms of momentum transport |25l 126] . The 
trace impurity transport dimensionless coefficients Ct,C u , and C p are properties of only the background 
turbulence, and can be calculated from 4 trace species (all boron in this case) with an (arbitrary) orthogonal 
set of gradients. Inserting the temperature and rotation gradients of the boron measured via charge exchange, 
we obtain a prediction for the steady state (r# = 0) density gradient which can be directly compared to the 
independently measured value. The predicted gradient is a combination of terms of different signs, requiring 
accurate treatment of symmetry breaking mechanisms and providing a sensitive test of the turbulence theory. 

The comparison of dimensionless parameters is particularly suited to quasilinear modelling, since only 
the shape and not the magnitude of the spectrum has an influence. The quasilinear results are obtained 
using a spectrum of 5 modes at kgpi — [0.2,0.3,0.4,0.5,0.6] in which the kgpi = 0.3 mode dominates (the 
choice of spectra is justified Sec. I4.3[) . Here, kgpi = kgp s ^jT e jTi where p s — c s /uj C i is evaluated using 
the field at the magnetic axis [j. The spectral shape is fixed the same for all cases and is compared with 
representative nonlinear simulations in Fig. [TU] (discussed further in Sec. 14.31) . This range of kg was used 
because wider ranges suffered from the problem of too many linearly stable modes. 

All the gyro-kinetic simulations include collisions (pitch angle scattering, with Z c g = l + (1.8el9/n e ) 6 ), 
full experimental flux surface geometry, and are performed in the frame of reference that rotates with the 
plasma, keeping both centrifugal and Coriolis forces (27] [28l E] and including the drive from the gradient in 
the plasma rotation [2H 1301 [31] ■ Electromagnetic fluctuations are neglected and are expected to give only 
minor impurity profile flattening for these cases [32] (which mostly have j3 e < 0.5%). 

The influence of neoclassical transport can be included by using the ion heat diffusivity Xi to 
independently normalize its contribution relative to the turbulent contribution [5]: 



The anomalous ion heat diffusivity xf n = Xi ~ Xi IS the difference between the power balance and 
neoclassical ion heat diffusivities computed using the transport code TRANSP [33], in which the neoclassical 
ion heat diffusivity \f C is calculated internally by NCLASS [31] ■ The ion heat diffusivity is chosen as the 
normalizing factor in preference to the total heat diffusivity because it avoids the need to calculate the 
radiation losses in the power balance analysis. 

| This definition of pi is inconsistent with our earlier definition of the thermal velocity, and is not the one actually used in 
GKW. In general geometry, the projection of kg onto the binormal coordinate uses the metric tensor [B], not all codes do this in 
the same way. 




(2) 




(3) 




(4) 



6 



10 u 



10" 



10" 



10" 



8 
s 



V 
V 




□ 

o 



V # 

o „ ^ y 



■ 


irn/GKi / ,,GK 
l RV B ' ' X i 







o 


_ RV Neo /x ar 


V 

+ 


D Neo /x an 



0.2 



0.4 0.6 
Mach B 



Figure 4. Relative importance of neoclassical and turbulent contributions to boron transport at r/a = 0.5, 
comparing each of the terms in Eq. [3] and also compared to the ratio of turbulent to neoclassical 
ion diffusivity. The quasilincar turbulent contributions are calculated for the kgpi = 0.3 mode using 
RVg K = D^ K (u'C u + (R/ Lt)Ct + C P ). Superscript abbreviations indicate the code used to calculate each 
quantity [gkw, TRANSP (an, inc. NCLASS), neo] . The four TEM points are also shown with an alternative 
normalizing factor (open quadrilaterals). 



The neoclassical RVg co is computed using the drift kinetic neoclassical code neo [9, 10, 11. Dimensional 
units are used to build the ratio RVg e ° /xt n , and as an independent cross- check, neo finds RVg eo /xf eo 
between and 4, with values consistent with the dimensionless comparison using the relation RV^ eo /xT = 
(RVg eo / xf e °) / (xf Ixf 1 )- The comparison of the neoclassical and turbulent contributions to the boron 
transport is shown in Fig. 01 Even though the neoclassical V/D may be quite large, the neoclassical 
contributions are at least five times smaller than the turbulent ones; when included in the predicted R/ L n B 
using Eq. [3]the differences are negligible, so for simplicity the neoclassical contribution is omitted in all other 
figures. Since the quasilincar gyrokinetic ion heat diffusivity may be under-predicted for TEM simulations 
(only four of the points), these points are also plotted with x2+l as the normalizing factor, but this does not 
alter the result that the turbulent contribution dominates. 

Input sensitivities of modelling 

Even when comparing dimensionless transport coefficients, the modelling of such a database requires very 
well diagnosed plasmas to generate inputs for the simulation. Since the mode properties largely depend 
on gradient quantities, very accurate profiles are required to determine the transport closely enough for a 
quantitative comparison. For each point in the database, the inputs to the modelling consist of plasma 
profiles in T e , Tj, n e , Vtor and the magnetic equilibrium including the q profile. The boron rotation is used 
for the main ion toroidal rotation (for all cases, the difference in deuterium Mach number is less than 0.03 
using Ref. [35]). Each profile fit must be carefully checked, so that bad calibrations and systematic features 
in the profile shape can be eliminated. Cross checking between independent diagnostics can also be useful; in 
this database Thomson scattering data was used where available to cross check the T e profiles from electron 
cyclotron emission (ECE) and the n e profiles from the IDA inversion. 

With this process completed, the most uncertain of the measured inputs to the simulations are the q 
profile (which is not well constrained by the magnetics measurements used to generate the equilibria) and 
the logarithmic density gradients (as discussed above). The temperature profiles are obtained from high 
resolution localized measurements (ECE and CX) of good quality, and do not suffer from the same level of 
uncertainty. 

To investigate the sensitivity of the modelling results to these uncertainties, a set of simulations was 
generated for the entire database, in which the most uncertain input parameters were varied within realistic 
estimates of their uncertainty (Fig. [SJ. The "nominal result" set is obtained using the full magnetic 
equilibria from the AUG function parameterization, while in order to conveniently vary the q profile inputs as 
independent parameters, the local geometry Miller parametrization l36j is also used (recently implemented in 
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GKW following the treatment in Ref. [37]; a benchmark with GS2 [35] is available [39 ]). Using this geometry, 
the q inputs were varied by ±0.5 (to a minimum of 1.1 for stability), and the magnetic shear s was varied 
by —0.3, +0.5 (which precludes any reversed shear configurations). In addition, the bulk plasma density 
gradients R/L ne were varied by ±0.8. 

The results of Fig. [S] confirm that the mode frequency at kgpi = 0.3 provides a good parameterisation 
of the particle transport over a large range of plasma conditions. We note that the main effect of varying 
the density gradient is to shift all the mode frequencies, with the transport coefficients lying along a single 
line. Varying the q profile, by contrast, alters both the frequency and the relationship between the mode 
frequency and the transport coefficients, increasing the predicted R/L n B most notably for ITG modes with 
a weak frequency, which were previously closer to the TEM transition. The largest enhancements to the 
convective C p pinch coefficient occur when the alterations to the q profile enhance the parallel compression 
part of the pinch (the term proportional to oc k 2 u>n in Eq. 10 of Ref. [22]), which for the ITG reinforces the 
E x B curvature pinch (oc |cj| 2 , first in Refs. [H32T]). The parallel compression pinch increases as magnetic 
shear increases and q decreases oc (s/qR) 2 ), and following Ref. [35] we refer to it as the "slab resonance". 
This provides the counterpart to the previously reported result of increased electron peaking with magnetic 
shear in the TEM regime [42] [43] , but we note that that effect is also partially due to the increase in the 
curvature drift frequency enhancing the inward thermodiffusion. These observations are consistent with the 
impurity transport mechanisms discussed in Ref. j22. 1441 [26] . 

The range of possible predicted R/L n B generated by varying the input q profiles and bulk species density 
gradients is used to generate an error bar on the 'nominal result' for each point (the result obtained using 
the default experimental profile fits (IDA for density and a function parametrization for the equilibrium). 
The sensitivities of the kgpi = 0.3 mode are used to assume the error for the quasilinear result including all 
5 modes (Fig. [5]). 

Given the sensitivity of the results to the slab resonance for the ITG mode, a further iteration of the 
sensitivity study is presented, using more realistic q profiles modelled with the ASTRA transport code coupled 
with the equilibrium solver spider 45.. Each shot is simulated in interpretive mode over a time window 
covering all the phases in the database. The current profile is allowed to evolve from an initial condition of 
(d/dr)(RE ( j > ) = 0, including torbeam [46 for EC current drive, and a simple sawtooth model to fix the 
profile inside the inversion radius. The resulting q profiles show good agreement with the measured inversion 
radius, and are the most accurate possible in the absence of a direct measurement. The comparison of these 
q profiles with those of the function parametrization in Fig. [TJindicates reduced q at mid radius, and reduced 
shear for the lower current cases (both conducive for the slab resonance). 

The simulation set with the improved q profiles has also been altered slightly to favour the ITG mode, 
relative to the nominal simulation parameters. For the 4 points which previously had TEM at kgpi = 0.3, 
Z e g is increased by 0.1, two points have R/L n decreased by 0.4, and one has R/Lxi increased by 0.4. In 
addition, the more physical friction and energy scattering terms in the collision operator have been added 
for the entire set. The resultant switch from TEM dominance at the largest scales is necessary for the 4 
previously TEM dominated cases to find sufficient levels of ion heat transport in comparison to the power 
balance calculation described in the next section. These slight changes to the inputs are enough to push 
the dominant modes at kepi — 0.3 to ITG, (some TEM modes persist at higher kg), to better match the 
heat flux ratio, and to allow the slab resonance in the impurity convection for these cases. The comparison 
between the simulations with both sets of q profiles and the measurement is shown in Fig. [5J 

4-2. Comparison of modelling and measurement, interpretation 

In this section, the results of the boron measurement and the quasilinear gyrokinetic modelling are compared 
and interpreted. In Figs. [6] and [8] both sets of results are plotted against the dominant experimental 
correlation of the boron Mach number, with the error bars generated from the sensitivity studies described 
in the previous sections included. Three main conclusions can be immediately drawn from these comparisons: 
First, both sets of modelling reproduce the experimental correlation, with the strength of the correlation 
sensitive to the q profile and the dominant mode. Second, nearly every point can be matched quantitatively 
within the estimated uncertainties, with the best agreement found for the improved q profiles and dominant 
ITG modes, which allow the most peaked cases to find the ITG slab resonance. Third, the modelling is 
only able to obtain hollow boron profiles in agreement with the measurement at high Mach number if the 
rotodiffusivc contribution to the boron transport is included. The quantitative agreement for the most hollow 
profiles is not perfect, an exact match would require either the rotodiffusive or thermodiffusivc component 
of the transport to be roughly a factor of two larger (see Sec. 14. 4p . 
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Figure 5. Sensitivity of heat transport and boron transport coefficients to input uncertainties in the safety 
factor and density profiles for the mode at kgpi = 0.3 and r/a = 0.5. From top: Ion/Electron heat flux 
ratio, Convective pinch, rotodiffusion, thermodiffusion, and predicted steady state R/L u b combining the 
coefficients using measured rotation and temperature gradients. Positive (negative) frequency u>n indicates 
propagation in the ion (electron) diamagnetic direction, signifying an ITG (TEM) mode. 

Figure 6. Comparison between predicted and measured boron gradients, using the nominal q profiles. 
(Top) Sensitivity of predicted R/L n g at r/a = 0.5 to the input uncertainties (for the dominant mode at 
kgPi = 0.3). The points and the legend correspond to those in Fig. \5\ (Mid) Comparison of the quasilinear 
(QL) predicted and measured values of R/L n g. The experimental errors bars are the same as those in Fig. 
[3] and the modeling error bars are generated from the sensitivity study of the top figure. The quasilinear 
predictions are calculated from a spectrum with its peak at kgpi = 0.3, which is compared with nonlinear 
(NL) simulations for selected points (discussed in Sec. I4.3D . (Lower) As mid, but excluding the contribution 
from rotodiffusion. 
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Figure 7. Comparison of example q profiles from AUG function parametrization and ASTRA. 

Figure 8. (Top) Quasilinear rotodiffusion coefficient for the improved q profiles from ASTRA with mild 
fTG mode biasing. (Mid) Comparison between predicted and measured boron gradients, for the nominal 
(circles) and improved (squares) quasilinear simulations. The measured points (diamonds) and the nominal 
simulations are the same as in Fig. [6] (Lower) Improved simulation case excluding the contribution from 
rotodiffusion, and comparison of the linear fits for each predicted case (the fit colors correspond to both the 
middle and lower panels). 
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Figure 9. (Upper) The mode frequency at kgpi = 0.3 is strongly related to the distance from the analytic 

Vl = L n /L T threshold for the ITG mode 07], where f? thrcs = 1.0 + y/l.O + (1 + T l /T e )(L n /_Rgfc fl p i ) 2 is 
derived in circular geometry. (Lower) Correlation between the mode frequency and boron Mach number is 
related to the AUG auxiliary heating systems (represented by the same symbol colors and shapes as in Fig. 
0. 



Given the demonstrated sensitivity of the modelling to the inputs, the combination of fluxes of opposite 
signs, and the sensitivity of the measurement to the deuteron profile, we find this agreement to be quite 
remarkable, building confidence that the local gyrokinetic model is correctly describing the dominant 
influences on the impurity transport. 

The success of the modelling allows us to understand the dominant experimental anti-correlation of 
boron peaking with the Mach number as a combination of correlations which connect the turbulence and 
the heating systems. The connections between the transport channels are framed in terms of the mode 
frequency, which we have demonstrated is the ubiquitous parameter for turbulent particle transport. 

The most fundamental relationship is that the auxiliary heating systems directly determine the 
temperature profiles, which set the characteristic frequency of the turbulence. For ITG modes this can 
be seen in the direct relationship between the mode frequency and rji = L n / Lt shown in Fig. [9] However, it 
can be more helpful to view this as a flux-driven relationship: The ratio of ion and electron turbulent heat 
fluxes Qi/Qe must match the power balance. This ratio has a direct relationship with the mode frequency, 
as demonstrated in Fig. El 

Since the turbulence dominates all transport channels, the selection of the mode frequency by the power 
balance has an immediate consequence for the electron density profile, in a relationship that is now well 
understood [48l [42l [26]. In the set of NBI heated, ITG dominated H-modes in the present dataset, this 
relationship manifests itself as an increased density peaking when electron heating is applied (as shown in 
Fig [I] and Ref. The density peaking increases as the electron heating shifts the turbulence from strong 

ITG towards the TEM transition (in the frequency domain) , which is the frequency region with the strongest 
convective particle pinch [321 122] ■ Of course, the turbulence is not determined by a single linear mode, so 
the transition from ITG to TEM is a gradual one, and feeds back on the heat flux relationship: Increased 
electron heating increases the prevalence of TEM even when ITG modes dominate; this increases the density 
peaking, which reduces rji and, therefore, the mode frequency further, a process that reinforces itself until 
the balance between ITG and TEM that matches the energy fluxes is reached. 

The rotation also has a correlation with the mode frequency (Fig. [5]), though the causes are less 
fundamental and more circumstantial: The NBI heating directly drives both rotation and ion temperature 
gradients (which affect the mode frequency through an increase of rji). It can also be seen in Fig. [9] that 
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this relationship depends on the beam geometry, since the most tangential NBI beam drives relatively more 
rotation for less increase in R/Lti at mid radius. In the more extreme geometry configuration of two 
oppositely directed NBI beams (as in DIII-D [3H]), this relationship could be broken altogether. Electron 
heating drives TEM, with the influence on the frequency and the electron peaking as discussed above (which 
also reduces rji). The increased density peaking has a consequence for the residual stress [33], which acts to 
reduce the rotation (Fig [1] and Ref. [!§]). 

Invoking these connections between the three channels of the bulk plasma, we are finally in a position to 
explain the anti-correlation between boron peaking and Mach number, as a combination of two effects: First, 
the (partially circumstantial) correlation between rotation and mode frequency (Fig. [5]) combines with the 
direct connection of the mode frequency to the impurity peaking (Fig. [5]) explaining part of the correlation. 
Second, the rotodiffusive contribution to the boron flux becomes important and at higher Mach numbers, 
because both the drive u' , and the C u coefficient increase (Fig. HI discussed further in Sec. I4.4[) . It is 
clear from Fig. [6l in which the rotodiffusive contribution C u u' is excluded, that both the correct turbulence 
frequency and the rotodiffusion are required together to explain the experimental observations. We also note 
that in the cases in which the relationship between frequency and rotation is weakest (due to the use of the 
more tangential beams), the rotodiffusive contribution is stronger, compensating the other terms to explain 
the near linear observed correlation. Since, as a trace species, boron has no back reaction on the turbulence, 
its transport can be considered an independent diagnostic of the turbulence theory. This is in contrast to 
the bulk transport channels in which the driving gradients must be necessarily consistent with the fluxes. 

In fact, the anti-correlation between boron peaking and Mach number is the most complex relationship 
in the database to describe, and the framework above also describes the other correlations visible in Fig. [2] 
The temperature gradients (and heating powers) have direct relationships with the mode frequency, which 
also connect to the density gradient. In addition, the boron density gradients are correlated with the electron 
density gradients because the same curvature pinch mechanism applies to all species. For impurities, the 
roto and thcrmodiffusive terms reduce the boron density peaking relative to the electron density peaking 
(for electrons, rotodiffusion is negligible, and thermodiffusion has opposite sign) 26 . 

4-3. Nonlinear results 

In this section we verify the quasilinear (QL) modelling against nonlinear (NL) simulations for three 
representative points from the database, at ub = 0.25, ug = 0.50 and ug = 0.74 (with the nominal q profiles) 
as shown in Fig. |5J As in the QL cases, each simulation contains four trace species to enable independent 
determination of each of the dimensionless boron transport coefficients. The nonlinear simulations are well 
resolved for energy fluxes (21 toroidal modes, 167 radial modes, 25 parallel grid points and 32 x 8 velocity 
grid points), but we have not excluded that the complete convergence in the impurity channels might require 
greater resolution. 

Each of the selected points have dominant ITG turbulence at larger scales, but the lower rotation cases 
have an increasing TEM character at smaller scales. It can be seen from the spectra that increased TEM 
activity shifts the spectral peak to smaller scales, and changes the shape of the spectra (Fig. [TU| . To enable 
robust simulation over a large database, our QL spectra is held fixed (more advanced techniques adapt the 
spectra depending on the linear growth rates |50j). 

For the two lower rotation cases, it can be seen in Fig. [S] that the differences between the NL and the 
QL results are no more significant than the input sensitivities of the simulation. The QL spectral peak is 
in good agreement with the NL one (Fig. [TU|) . and each of the impurity transport coefficients (Table [1} and 
predicted R/L n B match very closely. The heat flux ratios match less closely, because the QL spectrum used 
does not contain enough long wavelength modes which dominate the heat fluxes more. 

For the highest rotation ub = 0.74, the strongest ITG case, the spectral mismatch only partly explains 
the small discrepancy between the QL and NL results. In Table [IJ the QL result is re-weighted with the 
spectral peak at kgpi = 0.2 (dashed grey in Fig. fTO)) . which gives excellent agreement between C u and C p , 
but a remaining discrepancy in Ct of 0.07. Since the R/L n s result depends on this value multiplied by 
R/LiTh (= 5.8 for this case), accurate determination of the hollowest profiles requires an extremely precise 
determination of Ct- It can be seen that from the single mode values that Ct (and C u ) increase strongly 
at high kg for this case, so the exact spectral slope as well as the peak becomes important. 

As a cross-check on the simulations, the ion-electron heat flux ratios are also compared with the value 
from a power balance (PB) calculation in which the NBI ion/electron heating and equipartition power within 
r/a = 0.5 is computed by transp, and the radiation losses are assumed to be a factor Z e g/3.7 of the total 
auxiliary heating power within r/a = 0.5, and where all the ECRH power is centrally deposited (Sec. 14.11 
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Table 1. Comparison of dimensionless transport coefficients between nonlinear simulations, the quasilinear 
spectrum, and the quasilinear value for each individual mode. The three cases correspond to the three 
nonlinear points in Fig. [6] which use the nominal q profiles. Nonlinear values are time averaged for at least 
300 (R/v thi ) after saturation, and different averaging windows do not change any of the coefficients by more 
than 5% . All modes are ITG except where indicated TEM. 
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describes the Z e g approximation). Since this ratio is quite sensitive to the approximate radiation model, 
this quantity should only be treated as a qualitative check that the simulations are obtaining the relevant 
type of turbulence; indeed the values shown in Table [1] confirm that this is the case. 

For the four cases (for the nominal parameters) with dominant TEM even at the largest scales, nonlinear 
simulations find Qi/Q e much less than the power balance value (but similar to the QL result), which provides 
the motivation for the ITG biased set in Fig. [8] The nonlinear TEM runs are also much more sensitive to 
high-fcg energy pileup most evident in the vorticity spectrum |51 [ I52 j . and we have not yet found the correct 
dissipation to avoid this problem. Despite a saturated state being achieved with heat fluxes, the particle 
fluxes are particularly sensitive to these numerical difficulties, and give wildly unphysical results far from 
the QL results. The ITG cases shown include a small amount of k A hyper-dissipation which is sufficient to 
generate a clean vorticity spectrum (Fig. [TU1 defined here as n 2 e — hf). 

In summary, the quasilinear results using only the dominant mode at each scale, reproduce the nonlinear 
results for particle transport in ITG dominated turbulence down to a very precise level for every coefficient 
(when the correct spectral shape is used). In regions of transitions between turbulence regimes both the 
quasilinear and nonlinear simulations are less robust, making comparisons much more sensitive. 

4-4- Symmetry breaking and rotodiffusion 

The rotodiffusive contribution to the particle transport, C u u' , is driven by the gradient u' of the impurity 
rotation. The mechanism is a consequence of symmetry breaking mechanisms in the local gyrokinctic 
equation [25], such that the value of the coefficient C u depends on both the mode frequency (Fig. [5J, 
the rotation (Fig. [8]), and the nonzero fen developed by a mode due to the symmetry breaking (C' u oc 
2u — (Z Ab /2Ao){k^v t hi / oj d)) 0- Each symmetry breaking mechanism that can generate a momentum flux 
can also generate rotodiffusion. In our quasilinear simulations, three possible symmetry breaking mechanisms 
contribute to the rotodiffusion: The gradient in the bulk species rotation (responsible for diagonal momentum 
transport [29 J, the Coriolis force (which generates a momentum pinch [27] L and the up-down asymmetry of 
the flux surfaces [54] . The last of these is not present when the Miller geometry is used, and only has a very 
minimal effect at the r/a — 0.5 location studied. Symmetry breaking by E x B shear can also be considered 
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Figure 10. (Left) Comparison of nonlinear and quasilincar spectra of electrostatic potential, vorticity, and 
boron flux for species with no gradients for the u = 0.25 and u = 0.74 cases. (Right) Nonlinear particle 
flux for the trace boron species with {R/ Lx , Rj L n ,u'} = {0,0,1} for the u = 0.74 case, with and without 
E X B shear coupled for pure toroidal rotation [53| . For this trace, Vb = i~lbDb{Cu + C' p ). 



in the local model |53j , but as a coupling of many modes is only meaningful in nonlinear simulations. 

To assess the impact of this additional symmetry breaking mechanism, the two nonlinear cases at 
u b = 0.5 and ub — 0.74 were repeated including the E x B shearing that is required to produce a purely 
toroidal rotation gradient (in the cases without E x B shear, only the parallel projection of the rotation 
gradient is kept). The results, shown in Figs. 1151 and ITUl and Table[TJ demonstrate the Ex B shear can increase 
the rotodiffusion coefficient by up to 0.4 in the high rotation cases. This result demonstrates for the first 
time that each symmetry breaking mechanism that is relevant to momentum transport has its counterpart 
in the rotodiffusive particle transport, and that we have understood well the connections between these 
off-diagonal transport mechanisms [26j . Tests conducted by switching the sign of the plasma current in the 
simulation geometry indicate that for purely toroidal rotation, the E x B contribution to the rotodiffusion 
always acts to increase the rotodiffusion coefficient (for momentum transport this contribution always acts 
to reduce the effective momentum diffusivity [53]). 

In a global model, other symmetry breaking mechanisms can be present, and while formally smaller 
in terms of a = pi/L± small parameter expansion [25], are thought to be important in determining the 
non-negligible turbulent residual stress 43 . A leading candidate for relevant global symmetry breaking 
mechanisms is the profile shearing effect [S5] [5S1 [57] . The inclusion of these effects via global simulations 
might be expected to have a further influence on the rotodiffusive term, and is expected to increase the 
rotodiffusion coefficient in the ITG regime (in parallel to the inward directed contribution to the residual 
stress [43), which would improve the agreement for the most hollow R/L n s points at the highest Mach 
number. 



5. Conclusions 

In this work, high quality charge exchange measurements have been used to construct a database of boron 
density profiles for a range of H-mode plasmas in the ASDEX upgrade tokamak. In addition to better 
measurements, this database improves on previous work [5], by covering a larger range of auxiliary heating 
conditions, q profiles, and plasma rotation. This improved boron database allows quantitative comparison 
with modelling of the predicted logarithmic gradient of the boron density at mid radius, where neoclassical 
transport of boron is demonstrated to be negligible. The ability to perform quantitative comparisons 
represents a significant advance over the previous qualitative study, since it enables delicate conclusions 
to be drawn about the precise physical mechanisms at play. 

Quasilinear gyrokinetics is used to independently predict the turbulent contribution, by computing the 
ratio between turbulent diffusion and off-diagonal pinch, thermodiffusive and rotodiffusive components. To 
enable robust quantitative comparison, the effect of uncertainties in both the measurement and simulation 
inputs have been investigated. For the simulated dataset, the density profile and q profile have the largest 
input uncertainties, both of which influence the predicted boron gradient. Changes in the input density and 
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temperature gradients change the mode frequency and thus the particle transport, while changes in the q 
profile alter the contributions between parallel (slab) and perpendicular compression (curvature) resonances 
in the convective pinch for the ITG mode. These sensitivities highlight the need for excellent diagnostics and 
equilibrium reconstructions when making quantitative comparison with gyrokinetic modelling. Accurate q 
profiles are particularly vital for correct modelling of impurity transport channels due to the toroidal and 
slab resonances in the convective term. 

Within the investigated uncertainties, the modelling is able to quantitatively match the measured boron 
gradients, and reproduces the correlations observed in the experimental database. The primary correlation 
in the database, between boron density peaking and Mach number, is explained by the modelling as due 
to a combination of the turbulence regime and the rotation (resulting in rotodiffusion) , both of which are 
determined (and correlated) by the auxiliary heating power. 

The rotodiffusive contribution becomes important at the higher Mach numbers in the dataset, 
and is required to reproduce hollow profiles in agreement with the measured values. This constitutes 
an experimental confirmation of rotodiffusion, which is a consequence of the same symmetry breaking 
mechanisms which drive momentum transport. The inclusion of additional symmetry breaking mechanisms 
responsible for residual stress might increase the rotodiffusive contribution and further improve the 
agreement. Centrifugal effects were included in these simulations but do not significantly affect the results 
for the boron impurity at mid radius. 

In the absence of a direct measurement, the most accurate q profiles are obtained from an interpretive 
transport simulation coupled to an equilibrium solver. These q profiles allow the simulations to correctly 
determine the ITG parallel compression (slab) resonance in the impurity pinch, and give the best agreement 
between the boron predictions and measurement. The slab resonance is most significant for the more peaked 
boron profiles (the lower rotation cases), which are also nearest to the ITG / TEM transition. 

The strong correlation of all the boron transport coefficients with the turbulent mode frequency 
underlines its convenience as a ubiquitous parameter for characterizing turbulent particle transport. Since 
the ratio of ion and electron fluxes is also determined by the mode frequency and must match the power 
balance, simulations that match all transport channels consistently provide the most robust validation of the 
turbulent transport models. Indeed, in the context of uncertain inputs, iterating in the parameter space to 
match multiple transport channels could even be used to further constrain the measurement uncertainties. 
We believe the results presented in this work constitute a subtle and stringent test of the local gyrokinetic 
model, and build confidence in its predictive power with respect to the transport of light impurities. 
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